Load data

raw <- 
  bind_rows(
    read_excel(file.path("data", "20170727__data_export.xlsx")) 
  ) %>% 
  select(File, Nr., Start, Rt, End, # filename and retention times
         Ampl44=`Ampl 44`, Ampl45=`Ampl 45`, # amplitudes
         Area44=`rIntensity 44`, Area45=`rIntensity 45`, # areas
         Intensity=`Intensity All`, #peak intensities
         R13C = `R 13C/12C`, d13C = `13C/12C` # ratio and delta values
  )

Map peaks

### CHANGE MAPPING FILE NAME ###
data_matched <- map_peaks(raw, metadata_file = file.path("metadata", "smoky_hollow_peaks_July27.xlsx"), quiet = FALSE)
## INFO: 627 peaks in 26 files successfully assigned.
## WARNING: 95 peaks are unassigned or missing
## # A tibble: 95 × 22
##                                 File   Nr.    Start       Rt      End
##                                <chr> <dbl>    <dbl>    <dbl>    <dbl>
## 1  1387__170725_PTV_2uL_OG038 Ad.dxf    NA       NA       NA       NA
## 2  1387__170725_PTV_2uL_OG038 Ad.dxf    NA       NA       NA       NA
## 3  1387__170725_PTV_2uL_OG038 Ad.dxf    NA       NA       NA       NA
## 4  1387__170725_PTV_2uL_OG038 Ad.dxf    13 2019.149 2024.374 2030.435
## 5  1387__170725_PTV_2uL_OG038 Ad.dxf    18 2422.310 2426.281 2432.760
## 6  1387__170725_PTV_2uL_OG038 Ad.dxf    NA       NA       NA       NA
## 7  1387__170725_PTV_2uL_OG038 Ad.dxf    NA       NA       NA       NA
## 8  1388__170725_PTV_2uL_OG041 Ad.dxf    NA       NA       NA       NA
## 9  1388__170725_PTV_2uL_OG041 Ad.dxf    NA       NA       NA       NA
## 10 1388__170725_PTV_2uL_OG041 Ad.dxf    NA       NA       NA       NA
## # ... with 85 more rows, and 17 more variables: Ampl44 <dbl>,
## #   Ampl45 <dbl>, Area44 <dbl>, Area45 <dbl>, Intensity <dbl>, R13C <dbl>,
## #   d13C <dbl>, compound <chr>, is_ref_peak <chr>, map <chr>,
## #   Rt_target <dbl>, type <chr>, process <chr>, notes <chr>, depth <dbl>,
## #   n_matches <dbl>, n_overlapping <dbl>

Add standard values

standards <- read_excel(file.path("metadata", "gc_irms_indiana_A6.xlsx"))
data_w_stds <- data_matched %>% 
  filter(type == "standard", is_ref_peak == "no") %>% 
  left_join(standards, by = "compound") %>% 
  mutate(is_std = !is.na(true.d13C) | !is.na(true.d2H))

Evaluate standards

Visualize

data_w_stds %>% 
  ggplot() +
  aes(x = true.d13C, y = d13C, color = File) + 
  geom_smooth(method = "lm", se = FALSE, alpha = 0.5) +
  geom_point() +
  theme_bw() +
  theme(legend.position = "none")

Overview

std_corrections <- 
  data_w_stds %>% 
  group_by(File) %>% 
  do({
    m <- lm(d13C ~ true.d13C, data = .)
    data_frame(
      offset = coefficients(m)[1],
      slope = coefficients(m)[2],
      delta_ref_CO = -offset/slope,
      max_residual = max(summary(m)$residuals),
      rss = sum(summary(m)$residuals^2),
      r2 = summary(m)$r.squared
    )
  }) 
std_corrections %>% knitr::kable(d = 3)
File offset slope delta_ref_CO max_residual rss r2
1384__170725_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 10.543 0.997 -10.570 0.325 0.768 0.990
1386__170725_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 9.574 0.968 -9.895 0.434 2.247 0.969
1390__170725_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 7.805 0.900 -8.671 3.692 17.021 0.782
1391__170725_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 9.140 0.966 -9.457 0.645 6.846 0.912
1392__170725_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.106 1.043 -10.647 1.640 16.324 0.834
1393__170725_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 10.426 1.034 -10.083 2.153 24.046 0.771
1401__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 10.852 0.999 -10.858 0.213 0.387 0.995
1402__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 10.051 0.986 -10.198 0.674 6.453 0.919
1403__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 7.561 0.924 -8.178 1.320 18.578 0.777
1406__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 8.518 0.943 -9.030 0.839 9.965 0.871
1408__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 7.016 0.897 -7.818 1.051 8.840 0.873
1409__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.883 1.039 -11.437 0.613 2.863 0.966
1410__170726_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 6.952 0.899 -7.733 1.086 10.527 0.853
1418__170727_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 9.715 0.968 -10.033 0.420 2.041 0.972
1419__170727_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 8.094 0.935 -8.661 0.850 10.906 0.858
1420__170727_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 7.493 0.918 -8.161 1.123 13.791 0.822
1423__170727_A5_PTV_14ng@uL_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 7.747 0.938 -8.255 1.538 25.407 0.724

isotope values check by rentention time

data_w_stds %>% 
  ggplot() +
  aes(x = Rt, y = d13C, color = File) + 
  geom_point() +
  theme_bw() +
  theme(legend.position = "none") 

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

peak intensities check by rentention time

data_w_stds %>% 
  ggplot() +
  aes(x = Rt, y = Intensity, color = File) + 
  geom_point() +
  theme_bw() +
  theme(legend.position = "none") 

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

Data

message(sprintf("Infered reference tank isotopic composition: %.2f +/- %.2f",
                std_corrections$delta_ref_CO %>% mean(),
                std_corrections$delta_ref_CO %>% sd()))
## Infered reference tank isotopic composition: -9.39 +/- 1.17
# simplifed offset=ref gas correction (TODO: consider signal strengh (linearity correct), consider drift and evaluate standards between samples rather than average)
offset_avg <- std_corrections$offset %>% mean()
slope_avg <- std_corrections$slope %>% mean()

All

data_matched %>% 
  filter(type == "sample", is_ref_peak == "no") %>% 
  mutate(d13C.corrected = (d13C - offset_avg)/(slope_avg)) %>% 
  ggplot() +
  aes(compound, y = d13C.corrected, color = File) +
  geom_point() +
  theme(axis.text.x = element_text(angle = 90))

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

Against depth

data_matched %>% 
  filter(type == "sample", is_ref_peak == "no") %>% 
  mutate(d13C.corrected = (d13C - offset_avg)/(slope_avg)) %>% 
  ggplot() +
  aes(x = depth, y = d13C.corrected, color = File, size = Area44) +
  geom_point() +
  scale_x_reverse() + 
  facet_wrap(~compound, scales = "free") + 
  coord_flip()

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))